clear;
clear classes

%% Set parameters

invP = invParams;

%Mesh

cs = 25;
pad = round(cs*(1.2.^(1:15)))';
dx = [flipud(pad); ones(40,1)*cs; pad];
dy = [flipud(pad); ones(20,1)*cs; pad];
dz = [flipud(pad); ones(20,1)*cs; pad];

invP.mesh = meshDM(dx,dy,dz);
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations

xtx = sum(pad)+50:50:sum(dx)-sum(pad)-50;
ytx = ones(length(xtx),1)*(sum(dy)/2);
ztx = ones(length(xtx),1)*(sum(dz)/2)+30;

invP.loc = [xtx(:) ytx(:) ztx(:)];

%Time stepping
invP.disT = [250 1e-4];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%% Sigma and tau models

sigmaModel = zeros(invP.mesh.ncx,invP.mesh.ncy,invP.mesh.ncz) + 1e-8;
tauModel = zeros(invP.mesh.ncx,invP.mesh.ncy,invP.mesh.ncz);
etaModel = zeros(invP.mesh.ncx,invP.mesh.ncy,invP.mesh.ncz);

sigmaModel(:,:,1:invP.mesh.ncz/2) = 0.01;
tauModel(invP.mesh.ncx/2-1:invP.mesh.ncx/2+1,:,1:invP.mesh.ncz/2) = 0.1;
etaModel(invP.mesh.ncx/2-1:invP.mesh.ncx/2+1,:,1:invP.mesh.ncz/2) = 0.1;

invP.sigmaModel = sigmaModel(:);
invP.tauModel = tauModel(:);
etaModel = etaModel(:);

%%


[hs,~,~] = calcHs_TDIP(invP,etaModel,true);

d = projHs(invP,hs);

save d.mat d
